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Abstract 

A simple stochastic model which describes microtubule dynamics and explicitly takes into ac- 
count the relevant biochemical processes is presented. The model incorporates binding and un- 
binding of monomers and random phosphate release inside the polymer. It is shown that this 
theoretical approach provides a microscopic picture of the dynamic instability phenomena of mi- 
crotubules. The cap size, the concentration dependence of the catastrophe times and the delay 
before observing catastrophes following a dilution can be quantitatively predicted by this approach 
in a direct and simple way. Furthermore, the model can be solved analytically to a large extend, 
thus offering a valuable starting point for more refined studies of microtubules dynamics. 
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INTRODUCTION 



Microtubules (MT) are involved in key processes of cell functions such as mitosis, cell mor- 
phogenesis and motility. The building blocks of microtubules are a/3-tubulin heterodimers 
which can associate either laterally or longitudinally [l| . In biological systems, microtubules 
display unusual non-equilibrium dynamic behaviors, which are relevant for cell functioning. 
One such behavior, termed treadmilling involves a flux of subunits from one polymer end 



to the other, and is created by a difference of critical concentrations of the two ends 



a. 



In 



another behavior, termed dynamic instability, microtubules undergo alternating phases of 
elongation and rapid shortening [s!]. The two behaviors, treadmilling and dynamic instability 
result from an interplay between the polymerization and the GTP hydrolysis. 

The cap model provides a simple explanation for the dynamic instability: a growing 
microtubule is stabilized by a cap of unhydrolyzed units at its extremity, and when this 
cap is lost, the microtubule undergoes a sudden change to the shrinkage state, a so called 
catastrophe. The transitions between growth and shrinking can be described by a two state 



model with prescribed stochastic transitions 



5|. This model has lead to a number of 



theoretical and experimental studies {g-^, which have shown in particular the existence of 
a phase boundary between a bounded growth and an unbounded growth regimes. Although 
many features of microtubules dynamics can be captured in this way, this model remains 
phenomenological, because of the unknown dependence of the transition rates as function 
of external factors, such as tubulin concentration or temperature. 

To go beyond phenomenological models, one needs to account for the main chemical 
reactions occurring at the level of a single monomer [o], |lo|. These reactions can be assumed 
to occur between discrete states, and the corresponding transition rates can be observed 
experimentally. In this way, discrete models can be constructed, which capture remarkably 
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we., the mai,. dyuam.ea. features of siug.e act,,, or si„g.e ,„,c..otubu.e ma,„e„ts fl □ 
1141 ]. These discrete models have the additional advantage of being free from some of the 
limitations inherent to continuous models. 

The question of the precise mechanism of hydrolysis in microtubules or actin has been 
controversial for many years despite decades of experimental work. In the vectorial model, 
hydrolysis occurs only at the unique interface between units bound to GTP / ATP and units 
bound to GDP/AD P, while in the random model, hydrolysis can occur on any unhydrolyzed 



unit of the filament leading to a multiplicity of interfaces at a given time. Between these 
two limits, models with an arbitrary level of cooperativity in the hydrolysis have been 
considered (see for instance jl5, Q for actin and jsl for microtubules). The idea that the 
filament dynamics depends on the mechanism of hydrolysis in its interior or more generally 
on the internal structure of the filament has been recently emphasized and it has been given 
the name of structural plasticity As a practical recent illustration of that idea, the 

dynamical properties of microtubules can be tuned by incorporating in them GDP-tubulin 



in a controlled way 



In microtubules, many experimental facts point towards a mechanism of hydrolysis which 
is non-vectorial but random or cooperative. Studies of the statistics of catastrophes 
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already provided hints about this, but there are now more direct evidences. ^The observation 
of GTP-tubulin remnants inside a microtubule using a specific antibody [22] is probably one 
of the most compelling evidences. With the development of microfiuidic devices for bio- 
chemical applications, similar experiments probing the internal structure and the dynamics 
of single bio-filaments are becoming more and more accessible. Furthermore, it is now pos- 
sible to record the dynamics of microtubule plus-ends at nanometer resolution 23|, |2J], thus 
allowing essentially to detect the addition and departure of single tubulin dimers from micro- 
tubule ends. In view of all these recent developments, there is a clear need to organize all this 
information on microtubules dynamics with a theoretical model. Here, we propose a sim- 
ple one dimensional non-equilibrium model, accounting for the hydrolysis occurring within 
the filament. We show that this model successfully explains for many known experimental 
observations with microtubules such as: the cap size, the dependence of the catastrophe 
;ime versus monomer concentration and the delay before a catastrophe following a dilution 
-|21|. Our interpretation of this data confirms and goes beyond results obtained in a 
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recent numerical and theoretical study of the dynamic instability of MT 25 1. 

In vivo, the dynamics of microtubules is controlled by a variety of binding proteins, which 
typically modify the polymerization process. Here we focus on the physical principles which 
control the dynamic instability of microtubules in vitro in the absence of any microtubule 
associated proteins. Our model differs from previous attempts to address this problem, in 
that it is sufficiently simple to be analytically solvable to a large extend, while still capturing 
the main features of MT dynamics. 
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MODEL 



GTP hydrolysis is a two steps process: the first step, the GTP cleavage produces GDP-Pi 
and is rapid, while the second step, the release of the phosphate (Pi), leads to GDP-tubulin 
and is by comparison much slower. This suggests that many kinetic features of tubulin 
polymerization can be explained by a simplified model of hydrolysis, which takes into account 
only the second step of hydrolysis and treats tubulin subunits bound to GTP and tubulin 
subunits bound to GDP-Pi as a single specie 
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mm- 

This is the assumption which we 



make here. Therefore what we mean by random hydrolysis here is the random process of 
phosphate release, which as we argue, controls the dynamic instability of microtubules. 

Our second main assumption has to do with the neglect of the protofilament structure 
of microtubules. Protofilaments are likely to be strongly interacting and should experience 
mechanical stresses in the MT lattice. We agree that modeling these effects is important 
to provide a complete microscopic picture of the transition from the growing phase to the 
shrinking phase, since this transition should involve protofilament curling near the MT ends 

k 
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27l |. Here, we do not account for such effects, because as in Ref. [8|, we are interested in 
constructing a minimal dynamic model for microtubules, which would describe in a coarse- 
grained way the main aspects of the dynamics of this polymer. 

We also assume that the filament contains a single active end and is in contact with a 
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3: 



reservoir of subunits bound to GTP. The parameters of the model are as in Refs. 
the rate of addition of subunits U, the rate of loss of subunits bound to GTP, Wt, the rate 
of loss of subunits bound to GDP, Wn, and finally the rate of GTP hydrolysis r assumed 
to occur randomly on any unhydrolyzed subunits within the filament. In Fig. [H all these 
possible transitions have been depicted. We have assumed that all the rates are independent 
of the concentration of free GTP subunits C except for the on-rate lol], which is U = k^c. 
All the rates of this model have been determined precisely experimentally except for r. The 
values of these rates are given in table [11 

As a result of the random hydrolysis, a typical filament configuration contains many 
islands of unhydrolyzed subunits within the filament. The last island containing the terminal 
unit is called the cap. 
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FIG. 1. (a): Representation of the various elementary transitions considered in the model with 
their corresponding rates, U the on-rate of GTP-subunits, Wt the off-rate of GTP-subunits, Wd 
the off-rate of GDP-subunits and r the hydrolysis rate for each unhydrolyzed unit within the 
filament, (b) Pattern for a catastrophe with terminal units in the GDP state. 



TABLE I. Various rates used in the model and corresponding references. The conditions are that 
of a low ionic strength buffer. The value used for the rate of hydrolysis result from the analysis of 
the present paper. 



On-rate of T subunits at -|- end (^M ^) 3.2 [38] 



Off-rate of T subunits from + end 
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Off-rate of D subunits from -|- end 
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Hydrolysis rate (random model) 


r {s~^) 


0.2 





RESULTS AND DISCUSSION 



In this section, we obtain the nucleotide content of the filament within a mean-field 

Q, 14, 28|). We denote by i the 



approximation (for earlier references on this model, see 
position of a monomer within the filament, from the terminal unit at i = 1. For a given 
configuration, we introduce for each subunit i an occupation number r^, such that Xj = 1 if 
the subunit is bound to GTP and Tj = otherwise. In the reference frame associated with 
the end of the filament, the equations for the average occupation number are for i = 1, 



d{n) 

dt 



U{1 - in)) - Wt{ti{1 - T2)) + WD{T2il - Ti)) - r{n), 



(1) 
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and for i > 1 



+ Wd{{1 - ri)(r,+i - r,)) - r(r,). (2) 

In a mean-field approach, correlations are neglected, which means that for any i, j, (rjrj) is 
replaced by {Ti){Tj). At steady state, the left-hand sides of Eqs. [T]|2] are both zero, which 
leads to recursion relations for the (tj). Let us denote (ri) = q as the probability that the 
terminal unit is bound to GTP. The recursion relations have a solution of the form for i > 1, 

(3) 

where h = {U — q{WT + f))/{U — qWr)- Combining Eqs. [I]|31 one obtains q explicitly as 
function of all the rates as the solution of a cubic equation which is given in the appendix 
of Ref. 1^. The mean filament velocity (namely the average rate of change of the total 
filament length) is given by 

^ = (U -WTq-WD{l-q))d, (4) 

in terms of the monomer size d. At the critical concentration Cc, the mean velocity vanishes, 
which corresponds to the boundary between a phase of bounded growth for c < Cc and a 
phase of unbounded growth for c > Cc [l^ . The plot of this velocity versus concentration 
exhibits a kink shape near the critical concentration, which is not particularly sensitive to the 
mechanism of hydrolysis since it is present both in the vectorial and random model [ll|, |l^ . 



This kink is well known from studies with actin 



291] but has not been studied experimentally 



with microtubules except in Ref. [30l] in a specific medium containing glycerol. 

The distribution of the nucleotide along the filament length has a well defined steady-state 
in the tip reference frame at arbitrary value of the monomer concentration c. Using Eq. [3l it 
follows that (-Tj) = and therefore, the steady-state probability that the cap has exactly 

a length /, P^, is Pi = (IlLi (''"«) )(1 ~ This leads to the following expression: 

Pi = 6'a-i)/2^« (1 _ tfq-^ ^ (5) 
and, the corresponding average cap size is : 

(o = E^^' = E^'^'"^^v. (6) 

Z>1 Z>1 
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FIG. 2. Average cap size in number of subunits as function of the free tubulin concentration c in 
/iM. The hne is the mean-field analytical solution and the filled squares are simulation points. 



In figure [21 we show how this average cap size varies as function of the free tubulin 
concentration. The average cap becomes longer than approximatively one subunit above 
the critical concentration, Cc defined above, and which is about 7/xM for the parameters of 
table [H used here. At concentrations significantly larger than this value, the cap grows more 
slowly, as ^JiiU /2r as f/ — t- oo |l, 8|. In the range of concentration [0:100 //M], the cap stays 
smaller than about 47 subunits, which represents only 3.6 layers (or 28 nm). This estimate 
indicates that the cap is below optical resolution in the range of tubulin concentration 
generally used, which could explain the difficulty for observing it experimentally. 

A long standing view in the literature is that the cap could be as small as a single layer. 



as shown by experiments based on a chemical detection of the phosphate release [3l|- This 
view has been recently challenged by two experiments, in which the length fluctuations of 
microtubules were probed at the nanoscale, 23|,l2^. The interpretation of these experiments 



still generate debates [32 



33| . In any case, taken together these two experimental studies 



reported a highly variable MT plus-end growth behavior, which suggests that the cap size 
is a fluctuating quantity, larger than one layer but smaller than about 5 layers. We note 



that such a range is compatib 



from dilution experiments 



e with our prediction and agrees with the estimation obtained 
2l|. Furthermore, our stochastic model naturally incorporates a 



fluctuating cap size. Even if the cap is indeed below optical resolution, we note that t 



nis 



22|. 



does not rule out the possibility that it could be observed with the technique of Ref. 

In flgure O we also compare the predictions of the mean-fleld approximation with an 
exact simulation of the dynamics. We flnd that mean-fleld theory provides an excellent 
approximation of the exact solution when the free tubulin concentration is above the critical 
concentration, which corresponds to the conditions of most experiments 19|, |20| . Deviations 



can be seen between the exact solution and its mean-fleld approximation in flgure [2] but 
only below the critical concentration. Many other quantities of interest follow from the 
determination of the nucleotide content of a given sub unit, namely (rj), such as the length 
fluctuations of the fllament [l^l or the islands distribution of hydrolyzed or non-hydrolyzed 



subunits [1|, uM- These predictions should prove particularly useful in testing this model 
against experiments, since the island distribution of unhydrolyzed units or "remnants" will 
become accessible in future experiments similar to that of 
conditions. 



22 1 but carried out in in vitro 



Frequency of catastrophes and rescues vs. concentration 

One difficulty in bridging the gap between a model of the dynamic instability and ex- 
periments, lies in a proper definition of the event which is called a catastrophe, since the 
number of reported catastrophes is affected by several factors depending on the experimental 



23| 



conditions, such as for instance the experimental resolution of the observation 

Although a catastrophe manifests itself experimentally as an abrupt reduction of the total 
filament length, we choose to define it from the nucleotide content of the terminal region. 



Following closely Ref. [25|, we define a shrinking configuration as one in which the last 
units of the filament are all in the GDP state (irrespective of the state of the other units) as 
shown in figure [H The remaining configurations (with an unhydrolyzed cap of any size or 
when the number of hydrolyzed subunits at the end is less than A^) are assumed to belong 
to the growing phase. In such a two states description of the dynamics (with a growing 
and a shrinking phase), which is implicitly assumed in the analysis of most experiments, the 
catastrophe frequency fc{N) is the inverse of the average time spent in the growing phase, 



while the rescue frequency fr{N) is the inverse of the average time spent in the shrinking 
phase. It follows from this that the catastrophe frequency fc{N) can be obtained as the 
probability flux out of the growing state divided by the probabihty to be in the growing 
state. For instance for = 1, this flux condition is 

Ml)q = {WT + r)P,+rJ2Pj^ (7) 

i>2 

where the terms on the right proportional to Pi correspond to a transition of the terminal 
unit from the GTP to the GDP state, which can occur either through hydrolysis or depoly- 
merization of that unit, while the last term corresponds to hydrolysis of the terminal unit 
from cap states of length larger or equal than 2. We have derived the general expression 
of fc{N) in the case of an arbitrary as shown in Supporting Information (SI) Methods, 
and we have checked these results by comparing them with stochastic simulations using the 



Gillespie algorithm 



Ml- 



In the case of the vectorial model, the last term in Eq. [7] is absent and the catastrophe 
frequency is non-zero only below the critical concentration. The fact that catastrophes 
are observed in |19i] significantly above the critical concentration indicates that this data 
is incompatible with a vectorial mechanism. For this reason, we only discuss here the 
predictions of the random model. 

The catastrophe time Tc{N) = l/fc{N) is shown as function of growth velocity for N = 2 
in figure [3^, and as a function of the concentration of free subunits, c, for = 1 in figure 
|3]d. The growth velocity is simply proportional to the concentration of free subunits. For 
both plots, one sees that below the critical concentration which is in the range of 5-10 fiM, 
the catastrophe time is zero as expected since there is no stable filament in that region of 
concentration. Note that Tc{N) behaves linearly as function of c for = 2 but it behaves 



non-linearly for A^ = 1. Since the experimental data of 



19| shows a linear dependence, this 



comparison indicates that the data can be explained with the model for N = 2 but not for 



A^ = 1. The same observation has been made in Ref. 



25| . where the same data has been 



analyzed. Note however, in comparing this work with this reference the following differences: 



first, the model of Ref. 



25| neglects rescues and assumes that the duration of a catastrophe 



once started is zero while the present model includes rescues, and takes into account the 



finite rate of loss of GDP units. Secondly the results of Ref. 25| corresponds to the regime of 



high concentration of free subunits while the present model holds at any concentration even 
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FIG. 3. (a): catastrophe time Tc vs growth velocity for the N = 2 case, with the theoretical 
prediction (solid line) together with simulation points (a) and experimental data points taken 
from [19] for constrained growth (o) and free growth (•). (b): catastrophe time Tc vs free subunit 
concentration for the = 1 case, with the theoretical prediction (black solid line) and simulations 
(filled diamonds). In addition, the dash-dotted line and the dashed line represent respectively 
(T(l)) and (T(20)), which have been calculated using Eg. [271 



in the proximity or below the critical concentration. Thirdly the present approach leads 
to analytical results with the assumption that the filament has no protofilament structure 
while the results of Ref. j3] are numerical but that model includes a protofilament structure. 
Our analytical derivation of the catastrophe time confirms that the case = 1 differs in 
an essential way from the N > 2 case at high concentration. Indeed, the catastrophe time 
reaches a plateau when the concentration goes to infinity for = 1, while it goes to infinity 
for N > 2. This trend is already apparent in the figure [31 

In figure [21 we have used a value for the rate of hydrolysis r = 0.2, which is higher than 
that estimated in Ref. js] (there the estimate was 0.002). The reason is that the hydrolysis 
rate is a global factor which controls the amplitude of the catastrophe time, basically Tc{N) 
scales for an arbitrary N as The value r = 0.002 leads to a reasonable estimate 

for Tc{N) for = 1 (albeit with the wrong dependence on concentration), but if we take 
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seriously as we do here, the observation that only the definition with = 2 is compatible 
with the measured concentration dependence of the catastrophe time, then r must have a 
significantly larger value than expected, and 0.2 is the value that is needed for in order 
to match the experimental data. Finally, we also note that the scaling of Tc{N) as a power 
law of r means that large values of (such as > 2) can be excluded given the observed 
range of catastrophe times. 

We also show the distribution of catastrophe times calculated with the parameters given 
in table [H for = 1 and A^ = 2 in figure |H These distributions in both cases are essentially 
exponential (except at a very short time which is probably inaccessible in practice in the 



experiments), in agreement with the observations reported in Ref. 19|] with free filaments. 

One advantage of our microscopic model is that it can explain and predict different related 
aspects of the dynamic instability of microtubules. Specifically, it also allows to predict the 
statistics of rescue events when the polymer switches from the shrinking phase back into the 
growing phase. Assuming that the system reached a steady-state behavior, the frequency 
of rescues fr{N) can be calculated using fiux conditions similar to the ones used to obtain 
fc{N) (see SI Methods for more details). The corresponding expression is rather simple and 
it can be written as 

fr{N) = U + WDb''q. (8) 

We have carried out a complete numerical test of this frequency of rescues using stochastic 
simulations, which is shown in SI Fig. 1. 

Our model predicts that rescue events should be observable under typical cellular con- 
ditions and in experiments. However, surprisingly there is a very limited experimental 
information on rescues. The analysis of Eq. ([8]) might shed some light on this issue. At 
low concentrations of GTP monomers in the solution, when the rate U is small, the av- 
erage time before the rescue event, ~ 1/^5 might be very large. As a result, it might 
not be observable in experiments since the polymer with L monomers could collapse faster 
{Tcoiiapse — L/Wd) bcforc any rescue event could take place. At large f/, rescues are more 
frequent given that the polymer is in the shrinking state. But the frequency of the catastro- 



phes is very small under these conditions, the microtubule is 
phase. Therefore in these conditions, rescues are not observed 



most always in the growing 



3- 
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FIG. 4. (left)The distribution of catastrophe time (N=l) for different concentration values. C = 
9/iM(filled squares) and C = 12/iM(open circles), (right) The distribution of catastrophe time 
(N=2) for different concentration values C = 9//M(filled squares) and C = 12/xM(open circles). 
The distributions are normalized. 



First passage time of the cap and dilution experiments 

In dilution experiments, the concentration of free tubulin is abruptly reduced to a small 
value, resulting in catastrophes within seconds, independent of the initial concentration 
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35| . This observation is an evidence that the cap is short and independent of the initial 
concentration. The idea that the cap is short is also supported by the observation that 
cutting the end of a microtubule typically with a laser results in catastrophe. As we shall 
see below, all these well-known experimental facts about microtubules can be explained by 
the present model. 
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FIG. 5. Distributions of the first passage time of the cap for an initial cap of A; = 2 units, F2{t) as 
function of the time t, for various initial concentration of free monomers. The solid lines are the 
theoretical predictions deduced from Eq. (13) of the SI Methods after numerically inverting the 
Laplace transform, while the symbols are simulations. The circles correspond to a dilution into a 
medium with no free monomers, squares correspond similarly to a dilution into a medium with a 
concentration of free monomers of 2/xM, diamonds to 5/iM, and triangles to 9/iM. 

Here, we are interested in the time until the first catastrophe appears following the 
dilution. For simplicity, we take the definition of catastrophe introduced in the previous 
section for = 1, which means that a catastrophe starts as soon as the cap has disappeared 
(as shown in the previous section, one could extend this result to the more general case of 
an arbitrary A^). Let us then introduce Fk{t) the distribution of the first passage time Tk 
for an initial condition corresponding to a cap of length k, and a filament in contact with a 
medium of arbitrary concentration. As explained in SI Methods, it is possible to calculate 
analytically Fk{t), by a method recently used in the context of polymer translocation j^. 
After numerically inverting the Laplace transform of -Ffc(t), one obtains the distribution 
Fk{t) which is shown as solid lines in figure Ofor the particular case of A; = 2. As can be 
seen in this figure, the predicted distributions agree very well with the results obtained from 
the stochastic simulation in this case. 
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From the distribution Fk{t) we obtain its first moment, the mean first passage time of 
the cap {T{k)). As shown in SI Methods, we find that 

^ ^UWrJuKv) -UJn+i{y) 



where y = ^/WrJU, y = 2^UWt /r , n = {U + WT)/r, and the functions Jniy) are Bessel 
functions. The dependance of {T{k)) as a function of the initial size of the cap k is shown 
in figure El at small k, {T{k))_js essentially linear in k as would be expected at all k in the 



vectorial model of hydrolysis [13| , while here it saturates at large values of k (the value of this 
plateau can be calculated analytically but only for ?7 = see SI Methods). To understand 
this saturation, consider a cap which is initially infinitively large, then after a time of order 
1/r, the cap abruptly becomes of a finite much smaller size as a result of the hydrolysis of one 
unit at a random position within the filament. This feature will always happen irrespective 
of the monomer concentration, and indeed in figure El {T{k)) has a plateau for A; — )■ oo for 
all values of the monomer concentration. We note that such a behavior of (T(k)) as function 



of k has similarities with the case of non-compact exploration investigated in [37(], while the 
vectorial model of hydrolysis would correspond in the language of this reference to the case 
of compact exploration. 

Let us now turn to a practical use of this quantity for characterizing the dynamic in- 
stability. In the previous section, we calculated the catastrophe time Tc. We expect that 
this quantity is an average of {T{k)), and indeed we find for the case of = 1 that is 
bounded by (T(l)) and (T(20)) (the choice of 20 is purely illustrative) as shown in figure 
[3l The characteristic time observed in dilution experiments is another average of {T{k)). 
More precisely, let us denote {T{k))post as the first passage time in post-dilution conditions 
given that the initial length of the cap is k. The dilution time TdUuUon is then the average of 
{T{k))post with respect to the steady-state probability distribution of the initial conditions 
before the dilution occurs. In other words, 

Tdiiution = y^(T(fc))pogt-Pfc(predilution), (10) 

k 

where Pk{pre — dilution) is the stationary probability given in Eq. \5\m pre-dilution condi- 
tions. 

In the case that the final medium after dilution is very dilute, one can assume that 
the final free tubulin concentration is zero, which allows to simplify the general expression 
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FIG. 6. Mean first passage time {T{k)) as function of k for three values of the monomer concen- 
tration from bottom to top 0, 2 and 4/xM. The presence of steps in these curves is due to the fact 
that {T{k)) is only defined on integer values of k. Note the existence of a plateau for all values of 
the monomer concentration. 
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given in Eq. [27] as explained in SI Methods. Using Eq. [TOl one obtains the dilution time 
for the parameters of the table [T] which is shown in figure [71 The figure confirms that the 
dilution time can be as short as a fraction of seconds in this case. It is straightforward 
to extend this calculation to the case of an arbitrary value of the post-dilution medium 
{i.e for the case of a dilution of arbitrary strength) using the general expression derived 
in Eq. [271 As the amplitude of the dilution is reduced (by increasing the post-dilution 
concentration), the dilution time increases as well but the general sigmoidal shape remains, 
with in particular a plateau at concentrations above the critical concentration. The presence 
of these plateaux means that the dilution time is essentially independent of the concentration 
of the monomers in pre-dilutions conditions as observed experimentally. Note that the height 
of these plateaux scale with the hydrolysis rate. For instance, to explain the dilution times 



reported in [35|, one needs to use a smaller value of r as given in the table because of 
the use of the = 1 definition of catastrophe. Alternatively, just as in the calculation of 
the catastrophe frequencies, it is possible to keep the expected large value of r provided 
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FIG. 7. Dilution time (s) as function of free tubulin concentration (in ^uM) before dilution in the 
case that the post-dilution tubulin concentration is zero. Solid line is the mean- field prediction 
based on Eq. (TOland the symbols are simulation points. As found experimentally, the dilution time 
is essentially independent of the concentration of tubulin in the pre-dilution state, and the time to 
observe the first catastrophe is of the order of seconds or less. 

the N = 2 definition of catastrophe is chosen. Thus, complementary information can be 
obtained from the catastrophe frequencies and the dilution times. 

CONCLUSION 

In this work, we have explained several important features about microtubules dynamics 
using a model for the random release of phosphate within the filament. The results of 
our mean-field approach are analytical to a large extend. With this approach we could 
recover some well known features of MT dynamics such as the mean catastrophe time and 
its distribution or the delays following a dilution, but we have also investigated much less 
studied aspects concerning the cap size, the role of the definition of catastrophes (via the 
parameter A^) and the first passage time of the cap. The theoretical model and ideas 
presented in this paper for the case of microtubules could also apply to other biofilaments 

16 



such as actin or Par-M, for which the random hydrolysis model may be relevant as well. 
Furthermore, although the model describes a priori only single free filaments dynamics, it 
is also potentially useful for understanding constrained filaments, in the broader context of 
force generation and force regulation by ensembles of biofilaments. For this reason, it would 
be interesting to study extensions of the model to account for the various effects of MAPs 
on microtubules, which should shed light on the behavior of microtubules in more realistic 
biological conditions. We hope that this theoretical work will stimulate further experimental 
and theoretical studies of these questions. 
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Supporting Information 

DISTRIBUTION OF FIRST PASSAGE TIME OF THE CAP IN THE RANDOM 
MODEL 

Let us denote by Fk{t) the probability distribution of the first passage time of the GTP- 
tip (also called cap in the main text), for a cap which is initially of length k. This quantity 
obeys the following backward master equation, for /c > 1, 

dF (^~^ \ 

= U{Fk+i - Ffe) + WT{Fk-i -Fk)+r J] Fj - kF^j . (11) 

These equations are supplemented by the boundary condition Fo(i) = S{t). We will assume 
that the random walk followed by the cap is recurrent, which means here that the disap- 
pearance of the cap is certain, whatever the time it takes. That condition means that for 
all k>0, 

/•oo 

Fk{t)dt = 1. (12) 



^0 



'0 

We will make use of the Laplace transform of Fk{t) defined by 

/•oo 

Fk{s) = / e-''Fk{t)dt. (13) 
With this definition, the equations above take the following form, again for k>l 

k-l 

{s + WT + kr + U)Fk = UFk+i + WtF^-i + r ^ F„ (14) 

j=0 

with in addition the conditions Fo[s) = 1 and for all k > 0, Fk{s = 0) = 1, which follows 
from the normalization condition and the definition of the Laplace transform above. 

For the applications of this first passage time distribution to dilution experiments, we 
are interested mainly in calculating it using post-dilution conditions. In the case of a dilu- 
tion, the concentration of the free monomers following dilution is in general small. Let us 
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discuss separately the particular case where the concentration of the medium after dilution 
is zero in which case U = 0, and the general case of a dilution into a medium of prescribed 
concentration corresponding to [/ 7^ 0. 



Particular case of f7 = 

In this particular case, the recursion equations given in Eq. [H] are easy to solve. The 
solution is 

(fc-l m \ 
m=l j=l ^ ' / 

for > 1 with the convention that a sum over an index which ends at is void. The mean 
first passage time T(/c), is the first moment of Ffe(t) and thus satisfies {T{k)) = —dFk/dss=o. 
It follows that for A; > 1, 

In this particular case of [/ = 0, it is possible to derive an asymptotic form of this mean 
first passage time for — )• 00, namely (T) = limk^ooiTik)) ■ Indeed in this case, the sum 
can be written in terms of hypergeometric functions l|, |2l], and it reads 

(T) = ^-^F(1; — + 2;^). (17) 

Wt + r r r 

The expression of the mean-first passage time given in Eq. [16] can be used to obtain the 
delay before the appearance of the first catastrophe as explained in the main text. In this 
case, we find 

Td.^„^.o„ = ^Afc6'=(^-l)/2^^ (18) 
fe>i 

where 

1 ^ Wt 

~ Wt + t Wt+{3 + l)r ^^^^ 

This dilution time is shown in Fig. 7 of the main text. 



General case of ?7 / 

The solution to this general case is more involved but it can be obtained using Bessel 
functions (for a solution of a similar recursion see jsl). In a first step, we transform the 
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recursion of Eq. [H] using the difference variable Kk{s) = Fk{s) — Fk+i{s), which leads to 

{s + WT + {k + l)r + U)Kk{s) = UKk+i{s) + WtKu-i{s). (20) 



Then, we introduce the change of variable Kk{s) = y^gk{s) and we choose y = ^JWtJTJ in 
such a way that Eq. [20] takes the simpler form: 

s + f/ + VFT + (fc + l)r , , 
+ gk-As) = ^^^^ 

The solution to this equation can be obtained by comparing with the well-known identity 

Ju+i{x) + Ju-i{x) = —Mx), (22) 

X 

for Bessel functions. Thus, the solution has the form 

gki-s) = C J(^s+U+WT+{k+l)r)/r{y) , (23) 



where C is a constant and y = 2^UWt jr . The boundary condition given above for -Fo(s) 
leads to the following condition 

Uygi{s) + s 
''^'^ = s + Wr + r ^ 

which fixes the constant C. In the end, one obtains 

g^^g^ ^ sJis+U+WT+ik+l)r)/r{y) ^25) 

-UJ(^s+U+WT+r)/riy) + VUWTJ{s+U+WT)/riy) 

which satisfies in addition the required condition at s = namely that for all k > 0, 
gk{s = 0) = 0. With this expression, one obtains the Laplace transform of the first passage 
distribution of the cap, Fk{s) from 

k-l 

Fk{s) = l-J2y'9A^)- (26) 

j=0 

Although it is not immediately apparent, it can be checked that the particular case discussed 
above is indeed recovered by taking the limit f/ — t- of the general case. After using 
{T[k)) = —dFk/dss=o together with Eq. [221 one obtains the general expression for the mean 
first passage time of the cap {T{k)) given in the main text, which reads 

k—l 

VUWTJniy) - U Jn+i{y) 

where n = {U + WT)/r. 
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CATASTROPHES AND RESCUES FOR ARBITRARY N 



Catastrophes are associated with stochastic transitions between growing and shrinking 
dynamic phases. The microtubule is in the growing phase when it is found in one of poly- 
mer configurations with the unhydrolyzed cap of any size or when the number of already 
hydrolyzed monomers at the end is less than A^. We define Rk,i as a probability to be in the 
polymer configuration with I T monomers at the end that are preceded by D monomers 
(irrespective of the state of the other subunits), Qk,i as a probability to be in the polymer 
configuration with I D monomers at the end that are preceded by A; T monomers, and finally 
Sk^i is a probability that the last / monomers at the end are hydrolyzed except for the one 
subunit at position k counting from the end of the polymer. Formally these definitions can 
be also written as, 

Rkj = Prob{. ..D_^,T_^),Qk,i = Prob{. ..T_^,D_^), 

k I k I 

Sk,i = Proh{. . . p...T...D) . (28) 

I 

Note that the probability Pi to have the unhydrolyzed cap of exactly / monomers can be 
expressed as Pi = Ri^i, while the probability to be found in the growing phase is 

oo N-l 

Pgr = Y,Rl,l + Y,Q^^l- (29) 

1=1 1=1 

The simple mean-field theory assumes that the state of the monomer in the microtubule 
is independent of its neighbors, and it also estimates that the probability to find T or D 
monomer k sites away from the polymer is equal to h'^^^q or (1 — h^^^q) respectively, with 
the parameter h given by 

U-qWr ^ ' 

The probabilities defined in Eq. (128|) can be easily calculated yielding, 

Rk,i = 6'('-^)/^g''nV - Vq),Q,,i = 6W-i)/2^fc-Q(i _ 
j=i j=i 

fe-i I 

Sk,i = h'-\X{{l - V-'q) n (1 - ^''l)- (31) 
j=l j=k+l 

Then the probability to be found in the growing phase is 

N-l k 

P^^ = q + Y,b'qX{{l-V-'q). (32) 

k=i j=i 
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The frequency of catastrophes fdN) in steady-state conditions can be found from the 
fact that the total flux out of the growing phase, fcPgr, niust be equal to the flux to the 
shrinking phase, leading to the following equation, 

N 

fc{N)P,r = WtRn,i + r J2 Sk,N. (33) 

k=l 

Using Eqs. (13T!) and (!32|) . it can be shown that 

,,,n ^t^T nf.i(i - fr^o) + r zL, i'-' nf=i-»(i - ly-'i) W-zld - 
''^"^^ T^m^ItJT:^) ■ 

For = 1, we obtain a simple expression for the frequency of catastrophes, 

Ul) = WT{l-bq)+r, (35) 

while for = 2 it gives 

. .2) = W^t(1 - feg)(l - b\) +r[l-bq + b{l- q)] 

' 1 + 6(1 - g) ■ ^ ' 

A limiting behavior of the frequency of catastrophes for general A^ can be analyzed. For 
low concentrations of free GTP monomers in the solution, corresponding to m — 0, we have 
g — )• and 5 — )• 1 + r /wt, producing 

f,{N) ~ r + ^^-^ . (37) 

For large A^ and small hydrolysis rates [r /Wt ^ 1) the expression for the frequency of 
catastrophes is even simpler, 

UN) ^r + ^. (38) 

Another limit of interest corresponds to large concentrations {U ^ 1), where g — 1 and 
6 — i- 1, leading to fc{N) — )■ for all values of A^ > 2, while for A^ = 1 we have /c(l) r. 

This method of analyzing catastrophes can be also extended to calculating frequency of 
rescue events fr{N). The probability to flnd the microtubule in the shrinking phase is equal 
to 

7V-1 k 

Psh = l-Pgr = l-q-Yl b'qllil - V-'q). (39) 

k=i j=i 

The total flux out of this state is given by 



fr{N)P,h = UPsH + WdQi,n. (40) 
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which leads to the following equation 



fr{N) = U + 



(41) 



This expression can be further simplified to obtain the final result 



fr{N) = U + Wnb''q. 



(42) 



For all values of in the limit of f/ — )■ it yields — U, while for large U we have 



In addition, the average time before the catastrophe or before the rescue can be easily ob- 
tained by inverting the corresponding expressions for frequencies, namely, Tc{N) = l/fc{N) 



Numerical test of the predictions for rescues frequency 

We show here a comparison between the theoretical mean-field prediction for the rescue 
frequency given by Eq. |12] and results from stochastic simulations in figure O In the condi- 
tions of this figure, the filaments are sufficiently long and thus they do not collapse before 
rescue events occur. The theoretical mean-field predictions agree well with the simulations 
at concentrations of free monomers larger than the critical concentration. Deviations are ob- 
served at low concentrations near the critical concentration, in a way which has similarities 
with the deviations observed in Fig. 2. 
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and Tr{N) = l/fr{N). 
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FIG. 8. Rescue time (N=l) as function of concentration. The solid line is obtained from the 
mean-field theory. The data points (filled squares) are obtained from the simulations. 
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